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ABSTRACT 


One  strategy  for  discriminating  between  explosions  and  natural  events  depends  on  accurate  determinations 
of  event  locations,  including  focal  depths.  If  a  seismic  event  could  be  reliably  determined  to  have  a  focal 
depth  greater  than  a  few  kilometers,  one  could  be  confident  that  the  event  is  not  an  explosion.  But  to 
determine  focal  depths  accurately,  one  must  first  have  a  fairly  accurate  model  of  the  crustal  structure  in  the 
vicinity  of  the  event.  Unfortunately,  sufficiently  accurate  models  do  not  exist  for  many  regions  of  interest 
to  the  nuclear  explosion  monitoring  community.  Our  previous  work  focused  on  developing  and  evaluating 
strategies  for  locating  events  using  a  single  three-component  seismic  station;  velocity  models  were  obtained 
via  crustal  receiver  function  modeling,  and  waveform  correlation  methods  were  used  to  determine  focal 
depths  for  which  synthetics  fit  the  data  best.  However,  for  an  event  at  a  regional  distance  from  a  given 
station,  the  sampling  provided  by  the  teleseismic  phases  used  for  receiver  functions  is  not  ideal.  These 
waves  tend  to  approach  the  station  at  a  steep  angle,  sampling  just  a  narrow  cone  beneath  the  station. 

Better  sampling  is  provided  by  shear-coupled  PL  (SPL)  phases,  which  sample  the  crust  over  1000  km  or 
more  as  they  approach  the  station.  This  sampling  provides  a  better  lateral  average  of  the  crust  and  more 
closely  resembles  the  sampling  of  phases  emanating  from  seismic  events  at  regional  distances.  Our  current 
research  centers  on  modeling  SPL  phases  using  a  novel  modeling  algorithm  that  uses  the  reflectivity 
method  to  compute  synthetic  seismograms  while  holding  deeper  portions  of  the  mantle  fixed,  in  terms  of 
pre-computed  and  stored  reflectivity  and  transmission  matrices.  Layers  of  the  crust  and  upper  mantle  are 
allowed  to  vary  over  broad  ranges  and  the  entire  algorithm  is  powered  by  a  variant  of  simulated  annealing, 
a  global  optimization  method.  Using  the  results  of  all  trials  we  are  able  to  compute  the  posterior 
probability  distribution  for  model  parameters  and  thereby  assess  the  strength  of  the  constraints  placed  upon 
the  model  parameters  by  the  data. 

We  report  on  progress  in  developing  the  modeling  code  and  demonstrate,  through  synthetics  computed  for 
receiver  function  models  produced  for  China,  that  the  sampling  provided  by  SPL  is  significantly  different 
from  that  provided  by  receiver  functions.  Mismatches  between  synthetics  and  data  for  S-wave  windows 
including  the  Sp  converted  phase,  SsPmP,  and  SPL,  indicate  that  each  model,  while  it  may  represent 
structure  in  the  immediate  vicinity  of  a  station,  does  not  reflect  structure  over  a  broader  area  around  the 
station.  Given  the  sampling  of  the  SPL  phase,  trial-and -error  modeling  is  very  difficult  and  time- 
consuming.  These  results  confirm,  however,  that  Sp,  SsPmP  and  SPL  phases  carry  additional,  valuable 
information  that  can  help  constrain  crustal  and  upper  mantle  velocity  models. 

KEY  WORDS:  discrimination,  event  location,  crustal  structure,  wave  propagation 


OBJECTIVE 


Introduction 


“Receiver  function”  methods,  which  typically  use  the  reverberations  of  P  phases,  contained  in  the  P  coda, 
to  constrain  a  1-D  model  beneath  a  seismographic  station,  have  been  quite  successful  in  matching 
important  characteristics  of  existing  refraction  models.  Receiver  function  methods  allow  the  crust  beneath 
a  given  station  to  be  estimated  with  much  less  effort  and  cost  than  are  required  for  active-source 
experiments.  However,  because  receiver  function  methods  use  rays  that  arrive  at  the  station  with  steep 
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angles  of  incidence  (e.g.,  Zhao  and  Frohlich,  1996),  they  sample  only  a  narrow  cone  beneath  the  station. 

The  result  is  a  1-D  model  that  does  not  represent  a  broad  regional  average,  i.e.,  the  medium  through  which 
signals  from  a  later  event  of  unknown  provenance  and  location  would  propagate.  Furthermore,  the 
predominantly  vertically  propagating  reverberations  that  make  up  the  P  coda  reflect  steeply  off  the  crust’s 
internal  layers.  They  therefore  do  not  emulate  the  propagation  of  regional  phases,  which  reflect  at  more 
oblique  angles  or  are  refracted  by  these  layers.  This  is  a  problem  for  locating  and/or  determining  focal 
depths  of  small,  regional  seismic  events.  In  short,  because  of  the  data  they  use,  the  models  produced  by 
receiver  function  methods  may  be  inadequate  for  the  purposes  of  Comprehensive  Nuclear-Test-Ban  Treaty 
(CTBT)  monitoring. 

We  are  developing  a  new  method  for  modeling  the  crust  and  upper  mantle  that  retains  the  time  and  cost 
advantages  of  P-coda  receiver  function  methods  but  which  uses  types  of  data  that  are  more  appropriate  for 
CTBT  purposes:  Shear-coupled  PL  phases  (SPL),  Sp  phases  converted  at  the  Moho,  and  SsPmP.  SPL 
samples  the  crust  and  upper  mantle  in  the  vicinity  of  a  station  most  broadly  compared  to  Sp,  SsPmP  and  P. 
Modeled  simultaneously  (where  they  exist),  SPL,  Sp,  and  SsPmP  offer  the  potential  for  producing 
azimuthally  dependent  structural  models. 

What  is  shear-coupled  PL  P  Frazer  (1977)  and  Baag  and  Langston  (1985)  review  previous  work  on  SPL 
extensively;  we  will  summarize  aspects  that  are  relevant  for  our  proposed  CTBT  monitoring  research. 

The  conventional  S  phase  is  the  initial,  relatively  sharp  and  pulse-like  arrival  that  signals  the  beginning  of  a 
wavetrain  with  generally  longer  periods  and  normal  dispersion.  The  particle  motion  associated  with  the  S 
phase  is  rectilinear,  and  all  three  components  of  motion  are  in  phase.  The  dispersive  wavetrain  that  follows 
S  exhibits  prograde  elliptical  particle  motion  that  is  confined  to  the  vertical  plane.  Oliver  (1961)  named 
this  wavetrain  “shear-coupled  PL”  because  it  is  analogous  to  the  PL  wavetrain,  which  appears  between  P 
and  S  arrivals  at  regional  distances.  Oliver  (1961)  presented  a  theory,  based  on  the  observed  group  and 
phase  velocity  of  SPL  that  explained  the  phase  as  coupling  between  S  and  the  fundamental  leaking  mode  of 
Rayleigh  waves  in  the  crustal  waveguide.  According  to  this  theory,  shear  energy  radiated  by  an  earthquake 
(or  explosion)  travels  through  the  Earth  as  a  body  phase,  whereupon  it  impinges  upon  the  Moho. 

Afterward  it  travels  through  the  crustal  waveguide  as  trapped  P-waves  and  leaky  SV-waves.  The  only 
difference  between  a  PL  phase,  which  is  observed  at  regional  distances  from  a  source,  and  SPL  phases, 
which  are  observed  at  teleseismic  distances,  is  that  SPL  is  generated  by  a  shear  wave  impinging  upon  the 
Moho  at  regional  distances  from  the  observing  station.  Chander  et  al.  (1968),  Frazer  (1977),  and  Baag  and 
Langston  (1985)  essentially  validated  Oliver’s  hypothesis  by  presenting  computational  methods  for 
synthesizing  SPL.  However,  Chander  et  al.  (1968)  present  only  approximate  seismograms  and  do  not 
include  S  waves.  Frazer’s  method  does  not  include  the  effects  of  anelastic  attentuation,  nor  did  he  explore 
the  effects  of  near-source  structure  on  SPL  generation  and  propagation. 

While  previous  work  has  elucidated  the  origin  and  propagation  characteristics  of  SPL  (e.g.,  Oliver,  1961; 
Chander  et  al„  1968;  Alsop  and  Chander,  1968;  Poupinet  and  Wright,  1972;  Frazer,  1977;  Baag  and 
Langston,  1985)  attempts  to  use  a  portion  of  the  S  coda  to  model  the  cmst  explicitly  are  quite  limited  (e.g., 
Helmberger  and  Engen,  1974;  Swenson  et  al.,  1999).  Modeling  of  Sp  and  SsPmP  have  received  slightly 
more  attention  (e.g.,  Jordan  and  Frazer,  1975;  Langston,  1996).  This  disparity  is  understandable; 
synthesizing  SPL  is  computationally  intensive  (Frazer,  1977;  Baag  and  Langston,  1985).  Because  of  the 
difficulty  in  computing  waveforms,  the  inverse  problem  of  determining  structure  along  the  propagation 
path  from  observed  SPL  has  not  been  adequately  explored.  Our  objectives  are  (a)  to  evaluate  the 
usefulness  of  shear-coupled  PL,  Sp,  and  SsPmP  phases  for  modeling  crustal  and  upper  mantle  structure 
using  real  and  synthetic  data,  (b)  develop  a  waveform  inversion  technique  based  on  a  novel  implementation 
of  the  reflectivity  method  and  global  optimization  algorithms,  and  (c)  apply  the  inversion  technique  to 
model  the  crust  and  upper  mantle  beneath  China.  Figure  1  shows  examples  of  synthetic  seismograms 
computed  for  three  models  that  differ  in  their  upper  mantles  and  crusts.  One  can  clearly  see  the  sensitivity 
in  timing,  pulse  width,  and  relative  amplitudes  of  S,  Sp,  SsPmP,  and  SPL  phases.  Successful  modeling  of 
data  that  contain  these  phases  will  provide  strong  constraints  on  regional  P  and  S  velocity  models. 
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Figure  1 .  A  comparison  of  vertical-component  synthetic  seismograms  showing  S,  Sp,  SsPmP,  and  SPL 

waves  computed  using  the  models  shown  at  right.  These  seismograms  were  computed  using  our 
reflectivity  code  for  frequencies  between  0.001  and  1  Hz  and  a  dip-slip  dislocation  source  at  600 
km  depth. 
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RESEARCH  ACCOMPLISHED 


Development  of  the  Modeling  Technique 

The  reflectivity  calculation  involves  computation  of  reflectivity  matrices  for  a  stack  of  layers  as  a  function 
of  ray  parameter  (or  wavenumber)  and  frequency.  The  computation  of  reflectivity  responses  for  different 
ray  parameters  and  frequencies  is  completely  independent  of  each  other.  We  took  advantage  of  this 
independence  to  develop  a  reflectivity  code  that  can  be  run  on  parallel  computer  architectures.  We 
parallelized  our  code  over  the  ray  parameters,  i.e.,  to  each  node  we  assign  a  certain  number  of  ray 
parameters  to  compute.  The  number  of  ray  parameters  to  be  computed  can  be  determined  based  on  the 
processor  speed.  If  all  the  processors  are  of  the  same  speed,  they  are  assigned  equal  number  of  ray 
parameters.  At  the  end,  the  master  node  assembles  the  partial  responses  and  performs  the  inverse 
transformation  to  generate  synthetic  seismograms  at  the  required  azimuths  and  distances. 

We  used  MPI  for  message  passing  and  ran  our  code  on  a  PC  cluster  consisting  of  16  nodes;  each  node  is  a 
660-MHz  alpha  processor  with  8-MB  cache  and  1-GB  of  RAM.  A  Myrinet  interconnect  is  used  to 
communicate  between  nodes.  Figure  2  shows  timing  statistics  for  the  reflectivity  code  on  our  PC  cluster. 
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Figure  2.  The  statistics  of  an  example  run  to  compute  complete  seismograms  for  frequencies  from  0  to  0.5 
Hz  (number  of  frequencies  =  2000,  no  of  layers=192,  no  of  ray-parameters=6000)  on  the  DEC 
Alpha  cluster.  Note  the  near  linear  speedup  of  the  algorithm  with  the  increase  in  the  number  of 
processors. 


We  plan  to  obtain  further  increases  in  computational  speed  by  using  additional  processors  (we  have  been 
unable  to  use  all  16  processors  due  to  minor  hardware  difficulties),  tailoring  the  number  of  ray  parameters 
more  carefully  to  the  structure  being  modeled,  and  by  pre -computing  and  storing  reflectivity  matrices  for 
parts  of  the  mantle  that  will  not  vary  in  our  modeling.  The  high  frequencies  needed  to  compute  body 
waves  accurately  over  the  distances  to  which  SPL  is  observed  contribute  to  the  difficulty  in  computing 
SPL.  We  propose  to  essentially  remove  these  complications  by  storing  pre -computed  products  for  fixed 
parts  of  the  Earth  model. 

The  modeling  process  will  be  controlled  by  a  global  optimization  algorithm  called  Very  Fast  Simulated 
Annealing  (VFSA)  (e.g.,  Sen  and  Stoffa,  1995).  Simulated  annealing  (SA)  is  analogous  to  the  natural 
process  of  crystal  annealing  when  a  liquid  gradually  cools  to  a  solid  state.  The  SA  technique  starts  with  an 
initial  model  mo,  with  associated  error  or  energy  £(mo).  It  draws  a  new  model  mnew  from  a  flat  distribution 
of  models  within  the  predefined  limits.  The  associated  energy  £(m,,ew)  is  then  computed  and  compared 
against  E(m()).  If  the  energy  of  the  new  state  is  less  than  the  energy  of  the  initial  state,  the  new  model  is 
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accepted  and  it  replaces  the  initial  model.  However,  if  the  energy  of  the  new  state  is  higher  than  the  initial 
state,  mnew  is  accepted  with  the  probability  of  exp  ((ii(m;;ew )—  ^(mQ  ))/ T ),  where  T  is  a  control 

parameter  called  temperature.  This  rule  of  probabilistic  acceptance  (also  called  the  Metropolis  rule)  allows 
SA  to  escape  local  minima.  The  process  of  model  generation  and  acceptance  is  repeated  a  large  number  of 
times  with  the  annealing  temperature  gradually  decreasing  according  to  a  predefined  cooling  schedule.  A 
variant  of  SA,  called  Very  Fast  Simulated  Annealing  (VFSA)  speeds  up  the  annealing  process  by  drawing 
each  new  model  from  a  temperature  dependent  Cauchy-like  distribution  centered  on  the  current  model. 

This  change  with  respect  to  SA  has  two  fundamental  effects.  First,  it  allows  for  larger  sampling  of  the 
model  space  at  the  early  stages  of  the  inversion  (when  “temperature”  is  high),  and  much  narrower  sampling 
in  the  model  space  as  the  inversion  converges  and  the  temperature  decreases.  Second,  each  model 
parameter  can  have  its  own  cooling  schedule  and  model-space  sampling  scheme.  VFSA  therefore  allows 
for  individual  control  of  each  parameter  and  the  incorporation  of  a  priori  information. 


Uncertainty  Estimation:  PPD.  Posterior  Covariance  and  Correlation 

In  seismic  waveform  inversion,  more  than  one  model  can  often  explain  the  observed  data  equally  well  and 
trade-off  between  different  model  parameters  is  also  common.  It  is  therefore  important  not  only  to  find  a 
single,  best-fitting  solution  but  also  to  find  the  uncertainty  and  level  of  uniqueness  of  that  solution.  A 
convenient  way  to  address  this  issue  is  to  cast  the  inverse  problem  in  a  Bayesian  framework  (e.g., 

Tarantola,  1982;  Sen  and  Stoffa,  1995)  in  which  the  posterior  probability  density  function  (PPD)  is  the 
answer  to  the  inverse  problem.  The  PPD  is  proportional  to  the  product  of  a  likelihood  function  and  the 
prior.  It  can  be  shown  (e.g.,  Tarantola,  1982)  that  the  PPD  d(m  |d0^  )  is  given  by 

I  d obs  )oc  exP (-  E(m)) •  (6) 

Since  the  PPD  is  a  multidimensional  function,  one  generally  evaluates  marginal  PPD,  mean,  covariance 
and  some  higher  order  moments.  All  these  quantities  can  be  represented  by  the  following  multi¬ 
dimensional  integral 

/  =  Jt/m/(m)a(m|dofe).  (7) 

Thus,  we  are  now  left  with  the  problem  of  finding  efficient  methods  of  evaluating  a  multi-dimensional 
integral.  Importance  sampling  based  on  a  Gibbs  sampler  or  a  Metropolis  rule  (Sen  and  Stoffa,  1996)  can  be 
used  effectively  to  evaluate  these  integrals  and  to  estimate  PPD,  posterior  mean,  covariance  and  correlation 
matrices.  The  posterior  covariance  and  correlation  matrices  clearly  demonstrate  the  trade-off  between 
different  model  parameters.  Sen  and  Stoffa  (1995)  showed  that  multiple  VFSA  runs  with  different  random 
starting  models  could  be  used  to  sample  models  from  the  most  significant  parts  of  the  model  space.  These 
models,  when  used  to  evaluate  equation  (7)  result  in  estimates  that  are  fairly  close  to  the  values  obtained  by 
theoretically  correct  Gibbs  sampling.  Multiple  VFSA,  however,  is  computationally  very  efficient. 

Application  to  Data 

After  exploring  available  data  and  relevant  previous  studies,  we  decided  to  focus  first  on  China  (Figure  3), 
for  which  several  studies  of  crustal  receiver  functions  have  been  made.  The  fairly  thick  crust  beneath 
China  produces  good  separation  between  Sp,  S,  and  SsPmP  phases,  which  makes  China  an  ideal  location 
for  a  first  application  of  the  modeling  method.  We  selected  eight  deep  events  at  epicentral  distances 
ranging  from  22°  to  58°,  rotated,  filtered  and  windowed  the  data  around  the  S  arrival.  Figure  4  shows  data 
for  a  single  station  (BJT,  although  it  is  the  former  station  BJI  for  the  1 990  event),  including  the  results  of  a 
polarization  filter  (radial  component  multiplied  by  the  vertical  component),  which  aids  in  identifying 
phases.  Since  radial  and  vertical  motion  are  correlated  for  rectilinear  motion,  P-type  energy  shows  up  as 
positive  values  whereas  SV-type  energy  is  negative.  The  best  examples  of  Sp,  SV,  SsPmP,  and  SPL  are 
marked  with  arrows  in  Figure  4.  These  arrivals,  which  also  appear  clearly  in  our  sample  synthetics  (e.g., 
Figure  1),  are  highly  sensitive  to  the  structure  of  the  crust  and  upper  mantle  and  offer  strong  modeling 
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constraints  on  that  structure.  However,  due  to  the  nature  of  their  propagation,  modeling  these  phases  by 
trial-and-error  would  be  extremely  time-consuming  and  frustrating — if  it  is  feasible  at  all. 


Figure  3.  Map  of  China,  including  the  earthquakes  (stars)  and  stations  (triangles)  used  in  this  study.  See 
Epicentral  distances  range  from  22°  to  58°.  See  Table  1  for  earthquake  parameters. 


Table  1.  Events  used  in  this  study. 


Year 

Month 

Day 

Hour 

Min 

Sec 

Latitude  (°) 

Longitude  (°) 

Depth  (km) 

Mag 

1 

1990 

5 

24 

20 

9 

23.2 

-7.363 

120.363 

588.9 

6.4 

2 

1994 

9 

28 

16 

39 

51.7 

-5.786 

110.352 

637.5 

6.6 

3 

1994 

11 

15 

20 

18 

11.3 

-5.589 

110.186 

560.6 

6.5 

4 

1995 

8 

23 

7 

6 

2.8 

18.856 

145.218 

594.9 

7.1 

5 

1996 

5 

2 

13 

34 

29.0 

-4.548 

154.833 

500.0 

6.6 

6 

1996 

6 

17 

11 

22 

18.5 

-7.137 

122.589 

587.3 

7.9 

7 

1998 

2 

7 

1 

18 

59.5 

24.821 

141.746 

525.3 

6.4 

8 

1998 

8 

20 

6 

40 

55.8 

28.932 

139.329 

440.5 

7.1 
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Figure  4.  Vertical-,  radial-,  and  transverse-component  data  from  the  events  listed  in  Table  1  for  the  single 
station  BJT  (BJI  for  the  1990  event).  These  velocity  records  have  been  bandpass-filtered  with  a  5- 
pole  Butterworth  filter  with  corner  frequencies  at  0.01  and  0.5  Hz.  At  far  right  is  the  result  of  a 
simple  polarization  filter  consisting  of  the  vertical  component  multiplied  by  the  radial  component. 
Rectilinear  P-type  motion  is  indicated  by  positive  values  and  rectilinear  SV  motion  is  indicated  by 
negative  values,  which  aids  in  the  identification  of  phases. 

We  implemented  a  number  of  the  crustal  models  obtained  by  Mangino  et  al.  (1999)  and  computed  synthetic 
seismograms  via  reflectivity  for  deep  events  at  distances  of  30°-60°  to  compare  to  the  data.  These  P- 
velocity  models  were  converted  to  S  models  by  means  of  a  simple  proportionality  (1.8).  Synthetics  were 
computed  up  to  frequencies  of  0.5  Hz  and  were  convolved  with  the  Harvard  CMT  solutions  for  each  event. 
Our  synthetics  matched  some  aspects  of  the  data — generally  we  obtained  a  good  match  to  the  SV  wave  but 
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failed  to  match  the  Sp,  SsPmP,  and  SPL  phases  (Figure  5).  Even  this  level  of  agreement  is  impressive. 
Receiver  functions,  as  discussed  above,  model  reverberating  waves  that  sample  the  crust  differently  than  do 
the  shear  wave  converted  phases.  Second,  the  S  model  we  used  was  converted  simply  from  P  models. 
Lastly,  in  order  to  compute  synthetics  over  teleseismic  distances,  we  took  the  liberty  of  inserting  Mangino 
et  al.’s  models  into  PREM,  essentially  replacing  the  PREM  crust.  The  structure  of  the  uppermost  mantle, 
including  the  thickness  of  the  crust  and  the  velocity  contrast  across  the  Moho,  play  crucial  roles  in  the 
Sp/SV  time  separation,  waveshapes  and  amplitude  ratios.  Modeling  these  aspects  accurately  will  have  to 
await  the  completion  of  our  modeling  scheme. 


CONCLUSIONS  AND  RECOMMENDATIONS 

Although  we  believe  the  method  we  are  developing  holds  strong  promise,  we  have  several  concerns  to 
address.  Among  the  issues  that  must  be  settled  are  the  questions: 

Is  SPL  responsible  for  the  S  coda? 

An  important  assumption  underlying  this  research  is  that  at  teleseismic  distances,  SPL  waves  are 
responsible  for  a  significant  proportion  of  the  energy  in  the  S  coda.  Our  modeling  approach  would  be 
invalid  if  the  coda  structure  were  controlled  instead  by  fundamentally  different  processes,  such  as 
multipathing  or  scattering  from  lateral  heterogeneity.  However,  we  believe  that  the  latter  case  is  unlikely 
because  of  the  considerable  previous  success  of  investigators  such  as  Frazer  (1977)  and  Baag  and  Langston 
(1985),  who  generated  SPL  waves  to  reproduce  the  essential  features  of  signal  that  followed  the  S  phase. 
Moreover,  SPL  generally  contains  considerable  power  in  a  relatively  narrow  frequency  range,  and  even  in 
the  presence  of  scattered  energy  it  is  likely  that  usable  signals  will  emerge  with  suitable  filtering.  Finally, 
multipathing  may  be  responsible  for  coda  structure  at  some  stations.  However,  it  is  unlikely  that  this  is  true 
at  the  majority  of  stations;  identifying  the  near-station  geological  conditions  which  imply  that  SPL  does  not 
control  coda  structure  is  an  implicit  goal  of  this  project. 


Do  lateral  variations  in  structure  prohibit  1  -D  modeling  via  reflectivity? 

We  believe  that  the  relatively  long-period  nature  of  SPL  and  its  broad  averaging  properties  offer  an 
opportunity  to  obtain  an  accurate  1-D  average  of  even  laterally  heterogeneous  regions.  The  extent  to  which 
SPL’s  averaging  holds  for  strong  and/or  short -wavelength  lateral  variations  is  unclear,  but  it  is  likely  that 
large  regions  of  the  world  will  admit  1  -D  reflectivity  modeling  of  SPL  waves.  Our  method  would  provide 
an  accurate  (i.e.,  laterally -averaged)  1-D  model  of  the  propagation  path,  which  is  more  practical  for 
monitoring  purposes  than  a  3-D  model.  Our  contention,  which  we  will  evaluate  in  the  course  of  this 
research,  is  that  a  1-D  model  may  be  “good  enough”  to  locate  seismic  sources  in  much  of  the  Earth, 
provided  it  accurately  reflects  the  crustal  and  upper  mantle  structure  between  the  source  and  the  station,  not 
just  beneath  the  station. 
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Figure  5.  Comparison  of  data  from  the  1 1/15/94  Indonesian  event  (black  lines)  with  synthetics  (gray  lines) 
generated  using  Mangino  et  al.’s  (1999)  crustal  receiver  functions  inserted  on  top  of  a  PREM 
mantle.  Models  were  not  provided  for  ENH,  SSE,  and  XAN;  here  we  show  the  best-fitting 
synthetics  from  among  the  models  provided.  Note  that  the  SV  arrival  is  well  matched,  in  many 
cases,  but  the  Sp,  SsPmP,  and  SPL  phases  are  generally  poorly  matched.  The  good  SV  match  is 
impressive,  particularly  since  the  receiver  functions  consist  of  P  velocity  models — we  obtain  S 
velocities  through  a  constant  proportionality  (1.8).  The  poor  match  of  converted  phases  is  not 
surprising,  since  these  phases  sample  the  crust  more  broadly  than  do  the  P  phases  used  for  receiver 
functions  and  therefore  experience  different  structure.  This  comparison  illustrates  the  potential 
usefulness  of  these  phases  for  modeling  regional  structures  in  a  manner  most  useful  for  monitoring 
purposes. 


Can  we  avoid  significant  source-side  contributions  to  SPL? 

Baag  and  Langston  (1985)  elegantly  demonstrate  the  large  contribution  to  recorded  SPL  that  arises  from 
reverberations  near  the  source  for  shallow  earthquakes.  We  intend  to  model  SPL  propagation  in  the 
receiver  region,  which  means  we  need  to  isolate  receiver-side  from  source-side  contributions.  We  can  do 
this  most  effectively  by  restricting  our  applications  to  deep  events.  Consequently,  we  will  have  fewer 
events  at  our  disposal  and,  ultimately,  less  complete  coverage  of  the  world  than  if  we  were  to  use  both 
shallow  and  deep  events.  Applications  to  shallower  events,  in  order  obtain  more  complete  worldwide 
coverage,  can  be  explored  later. 
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Do  we  have  the  computing  power  to  invert  SPL  for  velocity  structure? 

The  process  of  inverting  SPL  to  find  crustal  structure  is  highly  nonlinear.  Even  though  we  will  use  state- 
of-the-art  forward  modeling  and  inverse  methods,  they  will  still  require  that  we  calculate  SPL  synthetics  for 
hundreds  or  thousands  of  crustal  models.  Figure  2  demonstrates  that  we  have  the  capacity  to  evaluate 
hundreds  of  models  in  a  day’s  time  and  that  computational  speed  increases  nearly  linearly  with  an  increase 
in  the  number  of  processors.  More  and  faster  processors  are  easily  obtainable,  so  computational  constraints 
are  not  a  problem. 

The  method  we  are  developing  will  be  especially  useful  for  CTBT  purposes  because  it  will  allow  us  to 
produce  structural  models  in  regions  that  (a)  have  little  or  no  seismicity,  (b)  average  over  a  broad  region  in 
the  vicinity  of  the  seismographic  station,  and  (c)  are  especially  well  suited  to  modeling  regional  phases,  due 
to  the  propagation  characteristics  of  these  converted  phases.  Models  with  these  characteristics  will  be 
useful  for  locating  small-magnitude  events  at  regional  distances  with  a  sparse  network  of  broadband 
stations,  such  as  the  International  Monitoring  System. 
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